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A  FORTRAN  III  (IBM  7090)  PROGRAM  FOR  THE  CALCULATION  OF  ELECTRON  WAVE 
FUNCTIONS  FOR  THE  SODIUM  ION 


This  report  describee  &  FORTRAN  III  (IBM  7090)  program  for  the  calcula¬ 
tion  of  sodium-ion  wave  functions*  This  has  been  done  in  view  of  future , 
theoretical  evaluations,  like  transition  probabilities  and  cross  sections, 
which  are  needed  for  guidance  and  interpretation  of  sod ium-aro- tunnel 
investigations* 

This  work  was  performed  under  NOL  Task  No.  FR-64* 
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PURPOSE* 


METHOD: 


A  FORTRAN  III  (IBM  7090)  PROGRAM  FOR 
THE  CALCULATION  OF  ELECTRON  WAVE  FUNCTIONS  FOR  THE  SODIUM  ION 


Determination  of  the  radial  part  of  the  ware  functions  for  the 
electrons  in  the  1-s,  2-s  and  2-p  shells  of  the  sodium  ion. 

The  Hartree-Fock  equations  are  used  for  this  in  a  simplified 
form,  so  that  one  deals  with  three  uncoupled  eigenvalue  equations. 
The  theory  and  the  application  of  the  results  will  be  discussed 
in  NOLTR  63-32. 


A.  Mathematical 


The  following  equations  are  solved  for  the  eigenvalues 

and  the  eigenfunctions  P(l,r),  P(2,r),  P(3,r)'  : 

2 

|^~2  +  f[y(r)+y0(l,l,r)]  -  e]L  j  P(l,r)  =  0 


I “2  +^yW+y  (2,2^)]  -  e 

l  dr  <• 


P(2,r)  =  0 


^2  +  |[y(r)+y0(3,3,r)  +  |  y2(3,3,r)]  - 


P(3,r)  =  0 


y(r)  =  11  -  2yo(l,l,r)  -  2yo(2,2,r)  -  6yo(3,3,r) 

OO 

v  k  /  k+1 

yk(n,m,r)  =  /(£)  P(n,p)P(m,p)dp  +  J  (?)  P(n,p)P(m,p)dp 

o  r 

with  the  boundary  conditions: 

r->0  P(l,r)^  r  P(2,r)^  r  P(3,r)^r2 

-e  .  r 

r  -►oo  P(n,r) _ .  e  n 

and  the  normalization  condition 

The  solutions  are  obtained  by  iteration:  Starting  with  a  first  guess 
for  ere2,e3  and  P(l,r),P(2,r)  ,P(3,r) 

l)  y(r)  and  yk(n,m,r)  are  calculated 


00 

j  P2(n,p)dp  =  1 
0 


where 


and 
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2)  Then  new  P(n,r)  are  calculated  by  integrating  the  differen¬ 

tial  equations  from  r  =  0  outwards.  This  is  repeated 
adjusting  until  the  absolute  value  of  P(n,r)  is  a 

minimum  at  an  arbitrarily  chosen,  fixed,  large  r-value.  The 
changes  in  are  effected  as  long  as  they  are  larger  than 

the  required  eigenvalue  accuracy  TE. 

3)  With  the  new  eigenvalues  Epe2,e3  and  tile  new 

eigenfunctions  p(l,r),  P(2,r)  P(3>r)  the  calculations  are 

started  again  in  l). 

oo 

This  is  repeated  until  j  (p0i^“*pnew) 2dr  "  TT>  the  prescribed 
vavef unction  accuracy.  0 
Numerical  Remarks: 


l)  Since  the  differential  equations  do  not  contain  the  first 

derivative,  Numerous  integration  procedure  is  used:  If  the  D.E. 

is  fM(r)  =  a(r)  •  f(r)  +  b(r)  one  can,  using  62f  .-(Ar)2[f 

J  J 

and  f,.-,  =  2fj  —  +  ,  derive 


J+l 


(Ar)- 


+  ^  [10 


!.Ma 

1  12  j+l 

2)  The  preceding  formula  shows  that  the  numerical  procedure  is 
simplified  if  the  values  of  the  independent  variable  r  are 
equidistant.  On  the  other  hand  the  accuracy  is  better  if  —  = 
constant.  Therefore  as  a  compromise  the  f-range  is  subdivided 
into  several  sets  and  in  each  of  these  sets,  the  values  are 
equidistant. 
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Since  the  eigenvalues  can  be  determined  numerically  only 
within  a  certain  accuracy,  the  pertaining  wave  functions  do 
not  decrease  exponentially  for  increasing  r,  but  their 
absolute  values  go  through  a  minimum  and  increase  then  exponen¬ 
tially*  Therefore  the  final  wave  functions  are  obtained  by 
integrating  the  differential  equations,  using  the  final  eigen¬ 
values  from  r  =  0  outward  and  from  large  i^-values  inward. 

The  outward  and  the  inward  solutions  are  then  matched  at  some 
point,  which  can  be  prescribed. 

Since  P(3*r)  enters  in  y(r)  with  the  coefficient  6,  P(l,r) 
and  P(2,r)  however  only  with  the  coefficient  r  ,  the  iteration 
procedure  is  speeded  up  by  determining  first  P(3>r). 

Programming 

Input:  in  the  following  order 
JF 

Format  (13)  ,  restriction  JF  k  20  .,  1  card 

JF  is  the  number  of  sets  into  which  the  range  of  the  independent 
variable  t  is  divided.  Each  of  these  sets  is  once  more  sub¬ 
divided  inbo  equally  spaced  intervals.  The  endpoints  of  these 
intervals  are  the  values  of  r  used  in  the  program  and  calculated 
by  the  program  in  the  subroutine  Radius. 

NF(J)  XF(J) 

Format  (14,  E14.6)  ,  restriction;  NF(JS(L))  ^  6  ,  J  =  1,...JF,JF  cards 
’NF(J)  is  the  number  of  intervals  into  which  the  Jth  set  is  sub¬ 
divided  . 

XF(J)  is  the  right  endpoint  of  the  Jth  set 
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JS(L) 

Format  (313)  ,  L  =  1, ...3  »  1  card 

JS(L)  is  the  number  of  the  r-set  in  which  the  backward  and 
the  forward  integrated  solution  for  the  Lth  wave  function 
will  be  matched. 

TE  TT  E(l)  E(2)  E(3) 

Format  (5E14.6)  1  card 

TE  is  permissible  absolute  error  for  eigenvalues. 

TT  is  permissible  sum  of  squared  deviations  of  wave  functions  . 
E(L)  is  the  (guess  of  the)  Lth  eigenvalue. 

P(L,N) 

Format  (3E14. 6)  ,  L  =  1, . . .3 •  M  =  1, ...NFF 
Restriction  NFF <  1000  ,  NFF  cards 

P(L,N)  is  the  Nth  value  of  the  guessed  lth  wave  function. 

NFF’  is  the  total  number  of  R-points  and  will  be  calculated 
in  Radius. 

DATE 

Format  (2A5)  1  card 

Output:  in  the  following  order 

Each  time  the  calculated  eigenvalue  E(L)  differs  from  the 
previous  value  of  E(L)  by  an  amount  smaller  than  TE  there  will 
be  printed: 
a)  P(L,N)  PB(N) 

Format  (2E18.9)  >  N  =  1,NFF 

P(L,N)  is  the  Nth  value  of  the  Lth  wave  function  obtained  by 
forward  integration. 

PR(N)  is  the  Nth  value  of  the  Lth  wave  function- obtained 
by  backward  integration. 
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b)  E(L) 

Format  (E18.9) 

E(L)  the  calculated  Lth  eigenvalue 
2)  At  the  end  of  the  complete  calculation  will  be  printed; 

a)  DATE 
Format  (2A5) 

b)  E(L) 

Format  (3E14.9)  ,  L  =  1,3 

E(L)  the  final  values  of  the  eigenvalues 

c)  TE,  TT 

Format  (2E14-6) 

TE  is  the  permissible  absolute  error  for  eigenvalues 
TT  is  the  permissible  sum  of  squared  deviation  of  wave 
functions, 

d)  R(N)  P(1,N)  P(2,N)  P(3,N)  Z(N)  ZNL(1,N)  ZNL(2,N) 

ZNL(3,N) 

Format  (8E14.6)  ,  N  -  1,NFF 

R(N)  is  the  Nth  value  of  the  independent  variable  r 
P(L,N)  is  the  Nth  value  of  the  final  Lth  wave  function 
Z(N)  is  the  Nth  value  of  the  sodium  ion  potential 
ZNL(L,N)  is  the  Nth  value  of  the  potential  for  the  Lth  wave 
function 

PROGRAM  FORM: 

A.  Type  of  Code  FORTRAN  7090 

B.  Type  of  Monitor  BELL 

C.  Number  of  Tape  Units  required  (other  than  System  Tapes);  None 
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CODING  INFORMATION: 

A.  Storage  Required 

1.  Regular:  3894  (Total) 

2.  Common:  19037  (Total) 

B.  Accuracy: 

TE  -  The  accuracy  to  which  the  eigenvalues  are  to  be  calculated 
can  be  prescribed, 

TT  -  The  sum  of  squared  deviations  of  iterated  wave  functions 
can  be  prescribed. 

C.  Timing: 

It  takes  about  7  minutes  to  calculate  70  values  for  each  of  the 

three  wave  functions,  when  the  prescribed  accuracy  of  the  eigen- 
-3 

values  is  TE  =  10  and  the  prescribed  sum  of  squared  deviations 
of  iterated  wave  functions  is  TT  =  lti’"’2  . 

D.  Additional  Subroutines  Required: 

This  program  uses  only  subroutines  on  the  library  tape  and  sub¬ 
routines  incorporated  in  this  program.  In  fact  this  program  has 
been  split  into  the  following  subroutines: 

1)  Main  Hartree  Storage  217  +  Common  is  the  Main  program. 

2)  Subroutine  BINT(L)  Storage  304  +  Common.  It  calculates  the 
Lth  wave  function  by  integrating  backwards  from  large  r  values 
to  small  ones  . 

3)  Subroutine  EIGEN(L)  Storage  283  +  Common.  It  determines  the 
Lth  eigenvalue. 

4)  Subroutine  FINT(L)  Storage  390  +  Common.  It  calculates  the 
Lth  wave  function  by  integrating  forward. 
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5)  Subroutine  OUTPUT  Storage  182  +  Common.  Prints  the 
final  results. 

6)  Subroutine  PRINT (L)  Storage  41  +  Common,  Prints  for  the 
I.th  wave  function  the  values  obtained  by  forward  integral 
tion  as  well  as  the  values  obtained  by  backward  integra¬ 
tion.  It  also  prints  the  Lth  eigenvalue. 

7)  Subroutine  RADIUS  Storage  59  +  Common.  Divides  the  r-range 
into  several  sets  and  covers  each  of  these  sets  by  equi¬ 
distant  points. 

8)  Subroutine  RINT(U)  Storage  137  +  Common.  Dimension  U(1000) 
calculates  the  integral  of  a  function  defined  on  the  points 
calculated  in  Radius.  U  are  the  values  of  the  integral, 
the  function  itself  is  transmitted  through  Common. 

9)  Subroutine  SSRL(L)  Storage  244  +  Common.  Calculates  certain 
integrals  of  the  Lth  wave  function 

10)  Subroutine  TEST  Storage  153  +  Common.  Tests  whether  the 
calculated  wave  functions  fulfill  the  required  accuracy. 

11)  Function  YB(A,D,L)  Storage  128  +  Common 

Dimension  A(10),  D(,2) 

Is  a  backward  interpolation  for  A,  to  get  from  one  set  of 
equidistant  r-values  into  the  preceding  one. 

12)  Function  YF(A,D,I)  Storage  81  +  Common. 

Dimension  A(7),  D(2) 

Is  a  forward  interpolation  for  A,  to  get  from  one  set  of 
equidistant  r-values  into  the  following  one. 
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13)  Subroutine  Y2(P2)  Storage  1088  +  Common 
Dimension  P2(l000) 

Calculates  some  integrals  of  the  wave  functions > 

14)  Subroutine  ZPI  Storage  101  +  Common.  Calculates  the 
potentials  for  the  wave  functions. 


USE; 

The  determination  of  these  wave  functions  was  undertaken  in  connection 
with  a  sodium-arc- tunnel  project.  With  their  aid  cross-sections, 
transition-probabilities,  line  strength  and  line  profiles  can  be  cal¬ 
culated.  This  in  turn  can  be  used  for  interpretation  of  experimental 
data. 

SPECIAL  FEATURES:  None 

RESTRICTIONS: 

A.  Error  Comments  and  Stops: 

The  program  calls  ENDJOB  if 

a)  NFF  >  1000  where  NFF  is  the  number  of  values  of  a  single  wave 
function. 

b)  NF(JS(L) )  <  6  where  NF(JS(I,))  is  the  number  of  equidistant 
points  in  the  r-range  in  which  the  forward  and  backward  in¬ 
tegrated  solution  of  the  Lth  wave  function  are  matched. 

B.  Range  of  Parameters: 

a)  NFF  <  1000  where  NFF  is  the  number  of  points  of  a  single  wave 
function. 

b)  JF  <  20  JF  is  the  number  of  sets  into  which  the  grange  is 
divided.  Each  of  these  ranges  is  again  subdivided  into  equal 


intervals. 
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APPENDIX  A 
Listing  of  Program 


FOR 

C  MAIN  HARTREE 

COMMON  JF*NF*XF»DX»NFF*TE*Tf*E*P*PF*R*S*RR*Z*ZNL»PB*JS*IS 
DIMENSION  NF(20)*XF(20)*DX(20)*E(3)*P(3*1000)*PF(3*1000)*R( 1000) * 
1S(3*1000)*RR(3*1000 ) »Z ( 1000) *ZNL ( 3 . 1 000 ) * PB ( 1000 > »JS (3 ) »I$<3  > 

10  READ  1000* JF 

20  READ  1001* (NF( J) »XF( J) »J«1.JF) 

21  READ  1003* { J S ( L ) *L3 1  ♦  3  ) 

22  DO  28  L  =  1  *  3 

23  JSSaJS(L) 

24  IF(6-NF(JSS) )25«25*180 

25  IS(L)»l+NF(JSS)/2 

26  JS1«JSS-1 

27  DO  28  J«1*JS1 

28  I S ( L ) ■ I S ( L ) +NF ( J ) 

30  NFF«1 

40  DO  50  J«1*JF 

50  NFF*NFF+NF< J) 

60  IF(NFF-IOOO)  70*70.180 

70  READ  1002*TE*TT*E< 1 ) » E  <  2 ) » E ( 3 ) ♦ ( (P(L*N) *1=1 .3) *N«1 »NFF) 

80  CALL  RADIUS 

90  DO  120  L«l*3 
100  CALL  SSRL(L) 

105  X*  SQRTF ( S ( L *NFF ) ) 

110  DO  120  N-l.NFF 
120  P(L*N)«P(L*N)/X 
130  CALL  TEST 
140  CALL  ZPI 
150  DO  160  L«l*3 
160  CALL  EIGEN (L ) 

170  GO  TO  90 

1000  FORMAT  ( 13 ) 

1001  FORMAT  (I4.E14.6) 

1002  FORMAT  ( 5E14.6/ (3E14.6 > ) 

1003  FORMAT  (313) 

180  CALL  SYSTEM 

END 
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FOR 

SUBROUTINE  BINT(L) 

COMMON  JF#NF*XF»DX*NFF*TE»TT*E*P*PF»R*S*RR»Z#ZNL*PB»JS*IS 
DIMENSION  NF ( 20  )  » XF ( 20 ) »DX ( 2 0 ) *E ( 3 ) »P ( 3 . 1000 ) »PF (  3 *  1000 ) . R ( 1000 ) # 
IS (3 » 10  00 ) *  RR ( 3*1000 >*2(1000) *ZNL(3*1000) *PB( 1000) *  JS ( 3 ) »IS(3) 

20  B«. 000000001 
30  W*DX< JF)**2/12. 

40  C«(ZNL(L.NFF-1)+E(L))*DX( JF)*#2 
50  X=C/2.+SORTF(C*(l.+C/2»>> 

60  PB(NFF)*B/(1.+X) 

70  PB(NFF-1)«B 
80  NA*NFF 
90  JA=JF-JS(L >+l 
100  DO  190  K*1*JA 
110  J»JF-*-l-K 
120  NF1-NF(J)-1 
130  DO  150  M*1 » NF1 
140  N*NA-1-M 

150  PB(N)»(PB(N+1 )*<2*  +  WMZNL(L»N+l)+E(L) 1*10. 

1  >+PB(N+2)*(W*(ZNL<L»N+2)+E< 

2L) >-l. ) )/<l.-W*(ZNUL»N>+E<L) ) ) 

160  NA»NA-NF(J) 

170  M*NA-1 

180  W«DX( J-l >**2/12. 

190  P8(M>«(PB(M+1 )*<2.+W*(ZNl(L*M+l)+E(U >*10. 

1  >+YF(PB(M+l) *DX( J-l ) »1)*(W* 

2 (YF(ZNL(L»M+1 )»DX( J-l) »3>+E(L) )-i. ) ) / ( 1 .-W* ( ZNL ( L *M ) ♦£ (U ) ) 

RETURN 

END 
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SUBROUTINE  EIGEN ( L ) 

COMMON  JF*NF*XF*DX*NFF*TE*TT»E*P*PF»R*S*RR»Z»ZNL*PB*JS*IS 
DIMENSION  NF(20)*XF(20)»DX(20)*E(3) »P 1 3*1000) *PF <  3  *  1000 ) *R < 1000) * 
IS ( 3.1000) *RR(3» 1000 ) *Z( 1000) *ZNL<3*1000) » PB ( 1000 ) ♦ JS ( 3 ) . I S ( 3 ) 

10  CALL  FINT(L) 

20  EL  =  E ( L  ) 

30  AL«P(L*NFF) 

40  E(L)*EL+«001 
50  CALL  FINT(L) 

60  IF(SIGNF(1»*AL)*P(L*NFF)  )  70.70*110 
70  EU*EL 
80  AU“AL 
90  GO  TO  210 

110  I F ( P ( L  * NFF  )  )  120*280*130 
120  I F ( AL-P ( L* NFF ) )  20*20*140 
130  1 F ( AL-P ( L*  NFF ) )  140  *20*20 
140  E ( L ) *EL-» 001 
150  GO  TO  50 
160  EO«E(L> 

170  E(L)»< AU*EL-AL*EU)/ (AU-AL) 

180  DE  »ABSF<EO-E(L) ) 

190  IF(DE-TE)  280.200*200 
200  CALL  FINT(L) 

210  I F ( P ( L . NFF ) )  220.280*250 
220  AL-P(L.NFF) 

230  EL*E< L ) 

240  GO  TO  160 
250  AU«P(L.NFF) 

260  EU«E(L) 

270  GO  TO  160 
280  CALL  B1NT(L) 

290  I SS« I S ( L ) 

300  X*P(L»ISS)/PB( ISS) 

310  DO  320  N»ISS*NFF 
320  P(L*N)»PB( N)*X 
330  CALL  SSRL(L) 

340  W«SQRTF(S(L.NFF) ) 

350  DO  360  N«1»NFF 
360  P(L.N)«P(L*N)/W 
370  CALL  ZPI 
380  CALL  PRINT ( L ) 

390  IF(ABSF(A-E(L) l-TE)  420*400*400 
400  A»E ( L ) 

410  GO  TO  10 
420  RETURN 
END 
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FOR 

SUBROUTINE  FINT(L) 

COMMON  JF*NF*XF»DX*NFF»TE*TT  * E *P »PF  ♦  R *S . RR ,Z *ZNL # PB * JS *  I S 
DIMENSION  N  F (20) »XF (20) »DX(20) *E(3) »P ( 3 » 1000 ) *PF  (  3  *  1000 ) »R( 1000)  * 
IS ( 3*1000) »RR(3*1000 ) *Z ( 1000 ) .ZNL ( 3 . 1000 ) . PB < 1000 ) * JS < 3  )  *  IS ( 3 ) 

10  FIRSTF(X)»X**(B+1.)*(1.-X  *1 1 . / < 1 .+B ) +X  **2 / ( 6 .+4.*B ) *< < 2 . -B I  * 

1121.-V)-11./(18.*(1.+B*9.) )*( <l.+B)*121.-(2.+5.*B)*V)*X  **3) 

20  B»FLOATF (L/3 ) 

30  W«=DX  ( 1 )  **2  / 12  • 

40  V=4.*( RR( 1 iNFF) /S ( 1 *NFF ) +RR ( 2 ,NFF ) /S< 2 .NFF ) +3 «*RR ( 3 ,NFF ) /S ( 3 *NFF ) 

1 )«2.*RR(L*NFF)/5<L.*NFF)-E(L> 

60  P ( L  » 1 ) »0 • 

80  P(l*2)»FIRSTF(R(2) ) 

100  P<L»3>*FIRSTF(R(3) > 

110  K*1 
120  KA»3 

130  DO  210  J“1  ♦  JF 
140  M=K+KA 
150  K«K+NF(J) 

160  DO  170  N«M»K 

170  P{L.N)«(P(L*N-1)*(2.+W*(ZNL(L»N-1)+E(U )*10. )♦  P ( L *N-2 ) * < W* < ZNL 
1 ( L ♦ N-2 ) +E ( L ) ) “1 • ) )/ (l.-W*(ZNL(L»N)+E<L>  )  ) 

180  IF(J-JF)  190*2 10 ♦ 190 
190  W=DX(J+1)**2/12. 

200  PIL.K+1  )«IP(L.K)*(2.+W*<ZNl(L»K)+E(L>  )*10.  )+YB(P(L.K-3  )  *DX(  J)  *3  )*( 
1W*(YB(ZNL(L»  K-3 )*DX(J) *3)+E(L)5-l.) ) / ( 1 .-W* ( ZNL ( L *K+1 ) +E ( L ) ) ) 

210  KA»2 
220  RETURN 
END 
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FOR 

SUBROUTINE  OUTPUT 

COMMON  JF*NF»XF*DX»NFF»TE*TT»E*P*PF*R»S»RR»Z*ZNL 

DIMENSION  NF(20)  *XF<20!  * DX  (  2.0)  »E(  3)  *P<  3*1000). PF(3  *1000)  *R<  1000)  » 
IS (3, 1000 ) *RR<  3.1000 ) »Z< 1000) #ZNL(3»1000) 

10  READ  1000 »  D1 #  D2 

20  PRINT1001*D1#D2,E(1 )  »E  ( 2 )  *E(  3  )  *TE»TT  »R(  1 )  *P  ( 1 » 1 )  ♦*>  (  2  *1 )  *P  (3. 1 ) » 

1Z< 1 ) »ZNL<1*1 ) »ZNL (2*1) »ZNL (3*1) 

30  M*  1 

40  DO  70  J«1*JF 
50  K«M+1 
60  M«M+NF(J) 

70  PRINT  1002  *  ( R  ( N  )  *  P  ( 1»N)*P(2*N)#P(3*N)»Z(N)»ZNL(1»N) *ZNL(2»N) »ZNL 

1  <3*N  )  *N«K»M) 

90  CALL  ENDJOB 

1000  FORMAT  (2A5) 

1001  FORMAT  ( 1H1 *36X *47HWAVE  FUNCTIONS  FOR  SODIUM  ION  BY  HARTREE  METHOD 
1//55X»2A5////17X»3HE1«E14*9,17X»3HE2=E14.9,17X»3HE3«E14*9/ 

2  25H  ACCURACY  OF  EIGENVALUE«E14.6*40X»26HACCURACY  OF  WA 
3VEFUNCTIONS»E14«6///3H  R13X *6HP IS ( R ) *8X . 6HP2S < R ) .8X  .6HP2P ( R ) »8X » 4 
4HZ( R ) ♦ 11X *6HZ 1 S(R  > #8X*6HZ2S< R ) *8X*6HZ2P(R)//8E14»6) 

1002  FORMAT  (1H  /(8E14«6)  ) 

1003  FORMAT  OE14.6) 

END 
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FOR 

SUBROUTINE  PRINTUI 

COMMON  JF*NF.XF»DX*NFF»TE»TT»E»P#PrtR*SiRR»Z*ZNLtPB*J5»IS 
DIMENSION  NF(20)»XF(20)*DX(20) »E(3) »P ( 3  *  1000 ) *PF ( 3 » 1000 ) »R<1000>* 
IS  (  3  *  1000)  »RR( 3*1000  )  » Z  (1000  )  *ZNU  3  » 1000  )  »PB  ( 1000  > » JS  ( 3  »  •  ISO  ) 

10  PRINT  1000  *(P(L»N)«PB(N)#N«1*NFF)»E(L) 

20  RETURN 

1000  FORMAT  (2E18*9) 

END 
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FOR 

SUBROUTINE  RADIUS 

COMMON  JF,NF,XF»DX»NFF»TE»TT,E,P|PF»R»S»RR»Z*ZNL 
DIMENSION  NF(20)*XF(20)»DX(20)»E(3)»P<3»1000)»PF(3»1000)  tRUOOO), 
15(3,1000) #RR( 3 i 1000) »Z( 1000) iZNU 3 *1000) 

10  K“1 
20  A»0» 

30  R ( 1 ) «0» 

40  DO  100  J»1,JF 
50  M»K+ 1 
60  <«K+NF(J) 

70  A«XF(J)-A 

80  DX( J)»A/FLOATF(NF(J)  ) 

A»XF(J) 

90  DO  100  N«M»K 
100  R(N)»R(N-1)+DX(J) 

110  RETURN 
END 
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FOR 

SUBROUTINE  R I  NT  <  U  ) 

COMMON  JF*NF*XF*DX*NFF*TE»TT»EtP»PF»R»S*RR*Z*ZiNL*PB»JS*lS 

1  »X1 

DIMENSION  NF(20) *XF(20)*DX(2Q) *E(3) *P ( 3* 1 000 ) * PF l 3  *  1000 > *R UOOO I » 
IS ( 3 *1000 )  *RR  (3*1000 )  *Z 1 1000)  #ZNL(  3*1000)  tPB ( 1000)  »  JSO )  •  ISO ) 

2  *U( 1000) *X1 ( 1000) 

A*DX (11/12# 

Xl< 1 >«0* 

U(1)»0* 

U(2)«A#(8.*Xl(2)-Xl<3> ) 

K*1 

DO  30  J«1»JF 
M»K+2 
K«K+NF(JI 
DO  20  N*M»K 
B«-Xl(N-2) 

C«8.*X1 (N-l ) 

D"5.#X1(N) 

20  U(N)»U(N-l)+A*IB«-OD) 

IF(J-JF)  25*30*35 
25  A«DX(J+1)/12* 

B»-YB(Xl(K-3  )»DX(J) *1) 

C»8.*X1(X) 

5*5,*X1(KM  ) 

U(K+1)*U(KM-A*(8-K>D> 

30  CONTINUE 
35  RETURN 
END 
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FOR 

SUBROUTINE  SSRLIL) 

COMMON  JF.NF*XF.DX.NFP*TE*TT.E*R.PF»R.S.RR»Z«ZNL*PB.JS»I$ 
DIMENSION  NFJ20LXF(20).DX<20>.E<3>*P(3*1000)»PF(3»1000>  .R(1000)» 
1S(3.1000).RR<  3.1000)  »Z<  1000)  »ZNU  3.1000)  »P0<  1000).  US<3).  ISO) 

10  R  ( 1 )  *  1 « 

30  S(L»1I«0« 

31  A«DXI1)/12. 

32  B»8.*P<L*2)**2 

33  C»-P(U3>**2 
40  S(L»2)*A*(8+C) 

90  RRU.l>«0. 

60  RR1L.2)  •  A*IB/R(2)-K7RO)  ) 

TO  K*1 

80  DO  1T0  J*1»JF 
90  M«K>2 
100  X«K4NF(J) 

110  DO  130  N*M.K  • 

111  B«-P(UN-2)**2 

112  C*8.»P(L.N-1)**2 

113  D«9.*P(L*N)**2 

120  5U.N>«SU»N-1>*A*(B+C«-D) 

130  RR(L*N)»RR(L»N-l)+A«(B/R{N-2)-K/R(N-l)+0/R(N) ) 

140  IFIJ-JF)  141. 170.190 

141  A*DX I J4-1 )  /12. 

142  B*-YB(P(U»K-3)*DX(J)»3H**2 

143  C*8.»Pa»tO»*2 

144  D«3.*P(L.K+1)**2 

150  5IUMl)»S<U»K)+A*<B+C-K)) 

160  RRI L»K*1 )*RR (L.K )+A*( B/R( K~1 )+C/R(K ) +D/R ( K+l ) l 
170  CONTINUE 
180  R ( 1 ) *0. 

190  RETURN 
ENO 
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FOR 

SUBROUTINE  TEST 

COMMON  JF#NF»XP#DX*NFF»T£#TT »E»P*PF»R»S»RR#Z#ZNL»PB#JSt  IS 
DIMENSION  NF120) »XF { 20 ) . DX ( 20  )  *E< 3 ) *P ( 3» 1000) *PF( 3  *  1000) »R( 1000 ) ♦ 
1S(3*1000)*RR(3*1000) *2(10001 »ZNU  3  *  1000  )  *PB  ( 1000 ) »JS(3)»IS(3) 

10  K«1 
20  T*0# 

30  DO  100  J»1*JF 
40  M*K+2 
30  K«K+NF(J) 

60  A»0# 

70  DO  80  N»M*K*2 

80  A»A+$Q  (N-2)*4#*$Q  (N-1I+SQ  I N  > 

90  A*A+Fl.OATF<NFU)-<NF<  J)/2)«2  )*<-$Q  (K-2)«'8**SQ  (K-1)+5#*SQ  (M) 

1/4# 

100  T»T+A*DX( J)/3# 

110  IF(T-TT)  120#120» 130 
120  CALL  OUTPUT 
130  DO  150  L»l#3 
140  DO  130  N“l #NFF 
130  PFU_#N)«(P(L»N)*PF(L*N))/2# 

160  RETURN 
END 
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FOR 

FUNCTION  VB<A»0»U 
DIMENSION  A ( 10 ) »D<  2  ) 

10  B»-DI2) 

20  Bl*  6  ♦  Dll) 

90  B2  •  Bl  ♦  Dll) 

40  89  •  82  ♦  D(l) 

90  C»6#*D(1)**3 

60  YB*fBl*B2*IA(3*l+l > *B3-A(1)*B  H-3#*B*83*{  A(L  +  l)*Bl-A(2*L+l  )*B2)  l/C 
70  RETURN 
END 


FOR 

FUNCTION  YF(A»D«l) 

DIMENSION  A (7) »D( 2 ) 

10  B«A(2*l+l)-2,*A<l+l)+A(l} 

20  CHI2*L+llH.»AIUll-3,**m 
30  X«D ( 1 ) /D ( 2 ) 

40  YF«(B*X+C>#X/2.+A(ll 
50  RETURN 
END 
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SUBROUTINE  Y2 ( P2  ) 

COMMON  JF*NF«XF»DX*NFF*TEsTT*E*P*PF*R*S,RR»Z»ZNL»PB»JS»IS 

ltXl 

DIMENSION  NF(20)»XF(20)*DX(20)*E(3) »P ( 3 *  1000 ) *PF ( 3 . 1000 ) »R< 1000) ♦ 
IS (3 *1000 1 *RR( 3*1000 ) »Z < 1000 ) »ZNL < 3 ♦ 1000 >♦ PB ( 1000 ) »JS(3) »  ISO ) 
2#P2llOOO) *X1 (1000) »U( 1000) 

DO  10  N«2*NFF 
10  Xl  (  N ) » (P ( 3  tN ) *R ( N  ) ) »*2 
CALL  RINT(U) 

DO  20  N»2»NFF 

20  X1(N)*P(3*N)**2/R(N )**3 
CALL  R I  NT ( P2 ) 

P2( 1 )*0« 

DO  30  N*2»NFF 
H*R  (  N  )  # *2 

30  P2(N)»U(N)/H+H*R(N)*(P2(NFF)-P2(N> ) 

RETURN 

END 


I 
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FOR 

SUBROUTINE  ZPI 

COMMON  JF*NF»XF*DX»NFF»TE»TT*E»P*PF*R»S»RR»Z*ZNL»PB»JS*IS 
DIMENSION  NF(20) »XF (20  I ♦  DX  I  20  ) *E( 3) »P(3ilOOO) *PF( 3 1 1 OOO ) »R( 1000) » 
IS (3 .1000) #RR(3*1000) iZ(1000) #ZNL  (  3  ♦  1000  )  *P0  <  1000 )  *  JS  ( 3  I  *  ISO) 

10  R ( 1 ) ■ 1  * 

CALL  Y2(Z> 

30  DO  70  N*2»NFF 

31  A»Y ( 1 »N ) 

32  B»Y(2»N> 

33  C*Y ( 3 »N ) 

40  EE  =-ll.*2.#(A+B+3.*C) 

41  D»2«/R(N) 

50  ZNL ( 1 *N ) *D*( EE  -A) 

60  ZNL ( 2  »N ) *D*( EE  -B ) 

70  ZNL ( 3 *N ) *D#( EE  ~C-*4*Z  (N)+1*/R(N)) 

80  R ( 1 ) *0 * 

90  RETURN 
END 
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MAIN  ROUTINE 


APPENDIX  B  -  FLOW  CHART 


4 
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BINT(L) 
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RADIUS 
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RINT(U) 
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Return 
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Y2(P2) 
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APPENDIX  C 
EXAMPLE: 


The  i^-range  was  chosen  from  0  to  10. 

It  was  divided  into  5  sets,  so  JF  =  5* 
These  sets  were 


from  0 

to  0.04 

0.01 

so  NF(1)  =  4  , 

XF(1) 

0.04 

0.2 

0.02 

NF(2)  =  8  , 

XF(2) 

0.2 

0.5 

0.05 

NF(3)  =  6  , 

XF(3) 

0.5 

1.2 

0.1 

NF(4)  =  7  , 

XF(4) 

1.2 

10 

0.2 

NF(5)  =  44, 

XF(5) 

the  forward  and  backward  integrated  solutions  get  matched  in  the  Interval 
0.5  <  r  <  1.2  so 


JS(l)  =  4  JS(2)  =  4  JS(3)  =  4 

The  permissible  absolute  error  for  the  eigenvalues  was  chosen  to  be  10  4 
so  TE  =  10-4  . 

The  permissible  sum  of  squared  deviations  of  wave  functions  was  chosen  to 
be  10~3 

so  TT  =  10~3  . 

The  guesses  for  the  3  eigenvalues  were  taken  as 

E(l)  =  81  E(2)  =  4-546948  E(3)  =  3.204463 

5 

Since  NFF  =  1  +  £  NF(J)  =  70  there  follows  three  columns  with  70  values 

j=l 

each  for  the  guesses  of  the  three  wave  functions. 

The  last  input  i3  the  DATE  01-25-1962 
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TRA 

5 

4  *04 

8  *2 

6  .5 

7  1.2 
44  10. 

4  4  4 


.1  E* 

“03 

.1  E-Q2 

81 

4.546948 

-0. 

-0. 

-0 

0 . 629350E 

00 

0 • 1 54058  E-00 

0 

400036  E-02 

0  •  11.2  86  3  E 

01 

0.274103E-00 

0 

140013E-0 1 

0.151784E 

01 

0.366138E-00 

0 

300027  E-0 1 

0 . 181301E 

01 

0.433163E-00 

0 

500045  E-Q 1 

0 . 219322E 

01 

0.5Q7191E  00 

0 

102009E-00 

0 . 23  573 1 E 

01 

0.517195E  00 

0 

16  3015  E-00 

0 • 237832E 

01 

0 .48 118  IE-00 

0 

231021 E-00 

0 . 230628E 

01 

0.412155E-00 

0 

302027  E-00 

0.217621E 

01 

0 .32012  IE-00 

0 

374034E-00 

0.201312E 

01 

0.2 1508  IE -00 

0 

444040E-00 

0 . 183502E 

01 

0  •  102038E-00 

0 

512046E  00 

0 . 165292E 

01 

-0. 150057E-0 1 

0 

57 8052 E  00 

0 . 12  2768  E 

01 

-0.30111 4E-QG 

0 

724065 E  00 

0 . 878488E 

00 

-O’. 5 552 09 E  00 

0 

842076E  00 

0.61 4341E 

00 

-0 •  763 288 E  00 

0 

9 3 3 0 8 4 E  00 

0 • 422235E 

-00 

-0 .92  3  348 E  00 

0 

99 7090 E  00 

0.287160E 

-00 

-0.103739E  01 

0 

103909E  01 

0 . 194108  E 

-00 

-0.111242E  01 

0 

1061 1 0 E  01 

0 . 880489E 

-01 

-0.116944E  01 

0 

1061 1 OE  01 

0 . 400222E 

-01 

-0.113943E  01 

0 

102109E  01 

0 . 180100E 

-01 

-0.106140E  01 

0 

956086E  00 

0 • 80Q445E 

-02 

-0.958361E  00 

0 

88Q079E  00 

0 • 300167E 

-02 

-0.846319E  00 

0 

79 9072 E  00 

0.100056E 

-02 

-0 • 736278 E  00 

0 

7 18065  E  00 

-0. 

-0 • 632238E  00 

0 

640058E  00 

-0. 

-0 • 455172  E-00 

0 

501045E  00 

-0. 

-0.319120E-00 

0 

38  5035  E-00 

-0. 

-0.220083E-00 

0 

292026  E-00 

-o.< 

-0 • 149056E-00 

0 

2 1 9020E-00 

-0. 

-0.101038E-00 

0 

163015  E-00 

-0  • 

-0.680257E-01 

0 

121011 E-00 

-0. 

-0.450170E-01 

0 

890080E-01 

-0. 

-0.300113E-01 

0 

650058E-01 

-0. 

-0 • 200075E-0 1 

0 

480043E-01 

-0. 

-0 . 130049E-01 

0 

35003 1 E-0 1 

-0. 

-0.900340E-02 

0 

260023  E-0 1 

-0. 

-0. 600226E-02 

0 

190017  E-0 1 

-0. 

-0.400151E-02 

0 

140013  E-0 1 

-0. 

-0.200075E-0 2 

0 

100009E-01 

-0. 

-0.914625E-03 

0 

69 148  9E— 02 

0. 

-0.120045E-02 

0 

488044E-02 

-0. 

-0 . 800305E-03 

0 

328030E-02 

0. 

-0.280106E-03 

0 

26  8024E-02 

0. 

-0.596273E-10 

0 

200016E-02 

-0  « 

0. 

0 

148Q13E-02 

0. 

0. 

0 

120011E-02 

-0. 

-0. 

0 

800072E-03 

0. 

0. 

0 

1 280025E-03 

0. 

0. 

0 

119220E-09 

0# 

0. 

0 
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0* 

•0# 

-0* 

0. 

Of 

Of 

Of 

0* 

Of 

0. 

-0. 

-Of 

0# 

0« 

Of 

0. 

-0. 

-Of 

Of 

0* 

Of 

0* 

0* 

Of 

Of 

-Of 

Of 

01-25-1963 
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CM  < 
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'O  u. 
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-O-OM-O 
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J-O'CO'COO' 
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UJ  UJ  U  UJ  UJ  UJ'  UJ  UJ  UJ  UJ  ID  UJ  UJ  UJ  UJ  U  UJ  UJ 

r-  O'  m  j  n  o  cm  -c  -o  o  r-  oo-r- jmo 
cMLA-ocovri'Of'i^J'N-cvji/i'—  -o  o  ao  o  -o  m 
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JW'O'-CMNJCMnll/lN'-inOCMO'O'M 

>jiniMucogin3,cnN'*-'UoinO'On 

mCMCMCM  —  —  —  —  —  —  —  —  —  —  OOCncD 


CO  o  o  o 


UJ  L1J  UJ  l-j 
JT  CM  —  CM 
■3-  O  O  O 

intnmn 
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